#### ===============================
# Title: Replication files for 'Analyzing International Organizations: How the Concepts We Use Affect the Answers We Get'
# Project: Code for figures 
# Author: Charles B. Roger; Sam S. Rowan
# Contact: croger@ibei.org; sam.rowan@concordia.ca
# Date: 1 April 2021
#### ===============================

library(tidyverse)
library(haven)

# Set working directory to 'Roger and Rowan RIO 2021' on your machine

#### Figure 1: Growth of IOs over time

## Load in data
informal_ios <- read_dta("data/input/informals_v3.dta")

## Create vector of years
panel <- data.frame(year = seq(1925, 2010, 1))

## Reshape informals data to get start and end years
informals_df <- informal_ios %>% 
  select(-state_name, -active) %>% 
  pivot_longer(names_to = "io", values_to = "member", acd:zc) %>% 
  mutate(member = ifelse(member<0, NA, member),
         member = ifelse(member>1, 1, member)) %>% 
  group_by(io, year) %>% 
  summarise(max_membership = max(member, na.rm=T)) %>% 
  mutate(max_membership = ifelse(is.infinite(max_membership), NA, max_membership)) %>% 
  filter(!is.na(max_membership)) %>% 
  group_by(io) %>% 
  mutate(n = row_number(),
         max = max(n)) %>% 
  filter(n == 1 | n == max) %>% 
  mutate(start_year = min(year),
         end_year = max(year)) %>% 
  distinct(io, start_year, end_year) %>% 
  arrange(start_year) %>% 
  mutate(n = 1) %>% 
  group_by(start_year) %>% 
  mutate(cumulative_informals_created = cumsum(n)) %>% 
  group_by(end_year) %>% 
  mutate(cumulative_informals_terminated = cumsum(n)) %>% 
  ungroup()

## Count cumulative number of informals up to a given year
informals_open <- informals_df %>% 
  select(start_year, cumulative_informals_created) %>% 
  group_by(start_year) %>% 
  filter(cumulative_informals_created == max(cumulative_informals_created)) %>% 
  rename(year = start_year) %>% 
  full_join(., panel, by = "year") %>% 
  arrange(year) %>% 
  mutate(informals_created = ifelse(is.na(cumulative_informals_created), 0, cumulative_informals_created)) %>% 
  ungroup() %>% 
  mutate(cumulative_ios = cumsum(informals_created)) 

## Count cumulative number of informals that have closed up to a given year
informals_closed <- informals_df %>% 
  select(end_year, cumulative_informals_terminated) %>% 
  group_by(end_year) %>% 
  filter(cumulative_informals_terminated == max(cumulative_informals_terminated)) %>% 
  rename(year = end_year) %>% 
  full_join(., panel, by = "year") %>% 
  arrange(year) %>% 
  mutate(informals_terminated = ifelse(is.na(cumulative_informals_terminated), 0, cumulative_informals_terminated)) %>% 
  mutate(informals_terminated = ifelse(year == 2010, NA, informals_terminated)) %>%
  ungroup() %>% 
  mutate(cumulative_ios_closed = cumsum(replace_na(informals_terminated, 0)))

## Merge and create an operating count
informals_join <- full_join(informals_open, informals_closed, by = "year") %>% 
  mutate(informals_operating = cumulative_ios - cumulative_ios_closed) %>% 
  select(year, informals_operating) 

## Load in formals data
formal_ios <- read_dta("data/input/formals_v3.dta")

## Reshape formals data to get start and end years
formals_df <- formal_ios %>% 
  pivot_longer(names_to = "io", values_to = "member", aaaid:wtouro) %>% 
  mutate(member = ifelse(member<0, NA, member),
         member = ifelse(member>1, 1, member)) %>% 
  filter(year<2011,
         year>1924) %>% 
  group_by(io, year) %>% 
  summarise(max_membership = max(member, na.rm=T)) %>% 
  mutate(max_membership = ifelse(is.infinite(max_membership), NA, max_membership)) %>% 
  filter(!is.na(max_membership)) %>% 
  group_by(io) %>% 
  mutate(n = row_number(),
         max = max(n)) %>% 
  filter(n == 1 | n == max) %>% 
  mutate(start_year = min(year),
         end_year = max(year)) %>% 
  distinct(io, start_year, end_year) %>% 
  arrange(start_year) %>% 
  mutate(n = 1) %>% 
  group_by(start_year) %>% 
  mutate(cumulative_formals_created = cumsum(n)) %>% 
  group_by(end_year) %>% 
  mutate(cumulative_formals_terminated = cumsum(n)) %>% 
  ungroup()

## Count cumulative number of formals up to a given year
formals_open <- formals_df %>% 
  select(start_year, cumulative_formals_created) %>% 
  group_by(start_year) %>% 
  filter(cumulative_formals_created == max(cumulative_formals_created)) %>% 
  rename(year = start_year) %>% 
  full_join(., panel, by = "year") %>% 
  arrange(year) %>% 
  mutate(formals_created = ifelse(is.na(cumulative_formals_created), 0, cumulative_formals_created)) %>% 
  ungroup() %>% 
  mutate(cumulative_ios = cumsum(formals_created)) 

## Count cumulative number of formals that have closed up to a given year
formals_closed <- formals_df %>% 
  select(end_year, cumulative_formals_terminated) %>% 
  group_by(end_year) %>% 
  filter(cumulative_formals_terminated == max(cumulative_formals_terminated)) %>% 
  rename(year = end_year) %>% 
  full_join(., panel, by = "year") %>% 
  arrange(year) %>% 
  mutate(formals_terminated = ifelse(is.na(cumulative_formals_terminated), 0, cumulative_formals_terminated)) %>% 
  mutate(formals_terminated = ifelse(year == 2010, NA, formals_terminated)) %>%
  ungroup() %>% 
  mutate(cumulative_ios_closed = cumsum(replace_na(formals_terminated, 0)))

## Merge and count the number of formal IOs operating
formals_join <- full_join(formals_open, formals_closed, by = "year") %>% 
  mutate(formals_operating = cumulative_ios - cumulative_ios_closed) %>% 
  select(year, formals_operating)

## Merge formal and informal counts, create total count
io_count <- full_join(formals_join, informals_join, by = "year") %>% 
  mutate(total_operating = formals_operating + informals_operating) %>%
  pivot_longer(names_to = "IO type", values_to = "Count", 
               formals_operating:total_operating) %>% 
  mutate(`IO type` = as.factor(`IO type`),
         Year = as.numeric(year),
         Count = as.numeric(Count)) %>% 
  mutate(`IO type` = case_when(`IO type` == "formals_operating" ~ "Formal IOs",
                               `IO type` == "informals_operating" ~ "Informal IOs",
                               `IO type` == "total_operating" ~ "Total IOs"))

## Plot
io_count %>% 
  filter(year>1959) %>% 
  ggplot(., aes(x = Year, y = Count, group = `IO type`)) +
  geom_line(aes(linetype = `IO type`, 
                color = `IO type`), 
            size = 1.8) +
  scale_linetype_manual(values = c("dotdash", "solid", "dashed")) +
  scale_color_manual(values = c("grey70", "black", "grey30")) +
  labs(x = "Year", y = "IOs operating") +
  geom_point(data = filter(io_count, year>1959 & `IO type` == "Informal IOs"),
             aes(Year, Count),
             size = 2.5, color = "black") +
  geom_point(data = filter(io_count, year>1959 & `IO type` == "Informal IOs"),
             aes(Year, Count),
             size = 1, color = "white") +
  theme(legend.position = "bottom", 
        panel.background = element_rect(fill=NA), 
        text=element_text(size=18),
        axis.text = element_text(size = rel(1.25)),
        panel.grid.major.y = element_line(colour = "grey80"),
        panel.grid.major.x = element_line(colour = "grey80"),
        panel.ontop = FALSE,
        axis.line = element_line(size = .5, colour = "black")) 
ggsave("data/output/figure1_iogrowth.pdf", units = "in", height = 6, width = 8)  

#### Figure 2: Scatterplots of membership in 4 years

# Define a common ggplot theme
scatter_theme <- theme(legend.position = "none", 
                       panel.background = element_rect(fill=NA), 
                       text=element_text(size=18),
                       axis.text = element_text(size = rel(1.5)),
                       panel.grid.major.y = element_line(colour = "grey80"),
                       panel.grid.major.x = element_line(colour = "grey80"),
                       panel.ontop = FALSE,
                       axis.line = element_line(size = .5, colour = "black"))

# Load in data 

io_membership <- read_dta("data/input/summaries_v3.dta")

# Make a plot for 1955
io_membership %>% 
  filter(year == 1955) %>% 
  ggplot(., aes(formals_sum_member, informals_sum_member)) +
  geom_point(size=5, shape=19, alpha=.85) + 
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  labs(title="1955", 
       x = "", y = "") +
  # x= "Formal IO memberships", 
  # y = "Informal IO memberships") + 
  scale_x_continuous(breaks = seq(0, 120, 30), limits = c(-0,125)) + 
  scale_y_continuous(breaks = seq(0, 80, 20), limits = c(-15,90)) +
  coord_cartesian(xlim=c(0,125), ylim=c(0,80)) +
  scatter_theme -> m_1955

# Make a plot for 1975
io_membership %>% 
  filter(year == 1975) %>% 
  ggplot(., aes(formals_sum_member, informals_sum_member)) +
  geom_point(size=5, shape=19, alpha=.85) + 
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  labs(title="1975", 
       x = "", y = "") +
  # x= "Formal IO memberships", 
  # y = "Informal IO memberships") + 
  scale_x_continuous(breaks = seq(0, 120, 30), limits = c(-0,125)) + 
  scale_y_continuous(breaks = seq(0, 80, 20), limits = c(-15,90)) +
  coord_cartesian(xlim=c(0,125), ylim=c(0,80)) +
  scatter_theme -> m_1975

# Make a plot for 1995
io_membership %>% 
  filter(year == 1995) %>% 
  ggplot(., aes(formals_sum_member, informals_sum_member)) +
  geom_point(size=5, shape=19, alpha=.85) + 
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  labs(title="1995", 
       x = "", y = "") +
  # x= "Formal IO memberships", 
  # y = "Informal IO memberships") + 
  scale_x_continuous(breaks = seq(0, 120, 30), limits = c(-0,125)) + 
  scale_y_continuous(breaks = seq(0, 80, 20), limits = c(-15,90)) +
  coord_cartesian(xlim=c(0,125), ylim=c(0,80)) +
  scatter_theme -> m_1995

# Make a plot for 2010
io_membership %>% 
  filter(year == 2010) %>% 
  ggplot(., aes(formals_sum_member, informals_sum_member)) +
  geom_point(size=5, shape=19, alpha=.85) + 
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  labs(title="2010", 
       x = "", y = "") +
       # x= "Formal IO memberships", 
       # y = "Informal IO memberships") + 
  scale_x_continuous(breaks = seq(0, 120, 30), limits = c(-0,125)) + 
  scale_y_continuous(breaks = seq(0, 80, 20), limits = c(-15,90)) +
  coord_cartesian(xlim=c(0,125), ylim=c(0,80)) +
  scatter_theme -> m_2010

# Paste plots together
scatter_fouryears <- gridExtra::grid.arrange(
  m_1955, m_1975, m_1995, m_2010,
  nrow = 2, ncol = 2,
  bottom=grid::textGrob("Formal IO memberships",
                        gp=grid::gpar(fontsize=30,font=8)),
  left=grid::textGrob("Informal IO memberships", rot = 90,
                      gp=grid::gpar(fontsize=30,font=8)))

ggsave("data/output/scatter_1955_2010.pdf",
       scatter_fouryears,
       width=16, height=12, units = "in")

#### Figure APP-IO context
io_context <- read_dta("data/input/context_v3.dta")

c_1981 <- io_context %>% 
  filter(year == 1981) %>% 
  ggplot(., aes(formal_io_context_v3, informal_io_context_v3)) +
  geom_abline(intercept = 0, slope = 1, color = "grey80") +
  geom_point(size=5, shape=19, alpha=.85) + 
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  labs(title="1981", 
       x = "", y = "") +
  # x= "Formal IO memberships", 
  # y = "Informal IO memberships") + 
  scale_x_continuous(breaks = seq(3, 7, 1), limits = c(3, 7.5)) +
  scale_y_continuous(breaks = seq(3, 7, 1), limits = c(3, 7.5)) +
  coord_cartesian(xlim=c(3,7.5), ylim=c(3,7.5)) +
  coord_fixed(ratio = 1) +
  scatter_theme

c_1985 <- io_context %>% 
  filter(year == 1985) %>% 
  ggplot(., aes(formal_io_context_v3, informal_io_context_v3)) +
  geom_abline(intercept = 0, slope = 1, color = "grey80") +
  geom_point(size=5, shape=19, alpha=.85) + 
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  labs(title="1985", 
       x = "", y = "") +
  # x= "Formal IO memberships", 
  # y = "Informal IO memberships") + 
  scale_x_continuous(breaks = seq(3, 7, 1), limits = c(3, 7.5)) +
  scale_y_continuous(breaks = seq(3, 7, 1), limits = c(3, 7.5)) +
  coord_cartesian(xlim=c(3,7.5), ylim=c(3,7.5)) +
  coord_fixed(ratio = 1) +
  scatter_theme

c_1995 <- io_context %>% 
  filter(year == 1995) %>% 
  ggplot(., aes(formal_io_context_v3, informal_io_context_v3)) +
  geom_abline(intercept = 0, slope = 1, color = "grey80") +
  geom_point(size=5, shape=19, alpha=.85) + 
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  labs(title="1995", 
       x = "", y = "") +
  # x= "Formal IO memberships", 
  # y = "Informal IO memberships") + 
  scale_x_continuous(breaks = seq(3, 7, 1), limits = c(3, 7.5)) +
  scale_y_continuous(breaks = seq(3, 7, 1), limits = c(3, 7.5)) +
  coord_cartesian(xlim=c(3,7.5), ylim=c(3,7.5)) +
  coord_fixed(ratio = 1) +
  scatter_theme

c_2005 <- io_context %>% 
  filter(year == 2005) %>% 
  ggplot(., aes(formal_io_context_v3, informal_io_context_v3)) +
  geom_abline(intercept = 0, slope = 1, color = "grey80") +
  geom_point(size=5, shape=19, alpha=.85) + 
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  labs(title="2005", 
       x = "", y = "") +
  # x= "Formal IO memberships", 
  # y = "Informal IO memberships") + 
  scale_x_continuous(breaks = seq(3, 7, 1), limits = c(3, 7.5)) +
  scale_y_continuous(breaks = seq(3, 7, 1), limits = c(3, 7.5)) +
  coord_cartesian(xlim=c(3,7.5), ylim=c(3,7.5)) +
  coord_fixed(ratio = 1) +
  scatter_theme

scatter_context <- gridExtra::grid.arrange(
  c_1981, c_1985, c_1995, c_2005,
  nrow = 2, ncol = 2,
  bottom=grid::textGrob("Formal IO context",
                        gp=grid::gpar(fontsize=30,font=8)),
  left=grid::textGrob("Informal IO context", rot = 90,
                      gp=grid::gpar(fontsize=30,font=8)))

ggsave("data/output/figure_app1_iocontext.pdf",
       scatter_context,
       width=12, height=12, units = "in")

#### Figure APP-IO joined

{
set.seed(20210513)
j_1975 <- io_membership %>%
  filter(year == 1975) %>% 
  ggplot(., aes(formals_sum_joined, informals_sum_joined)) +
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  geom_abline(intercept = 0, slope = 1, color = "grey80") +
  geom_jitter(width = 0.25, height = 0.25, size = 4, alpha = 0.65, shape = 19)+
  labs(title="1975", 
       x = "", y = "") +
  scale_x_continuous(breaks = seq(0, 10, 2), limits = c(0, 10)) +
  scale_y_continuous(breaks = seq(0, 8, 2), limits = c(0, 8)) +
  coord_cartesian(xlim=c(-1,11), ylim=c(-1,11)) +
  coord_fixed(ratio = 1) +
  scatter_theme

set.seed(20210513)
j_1985 <- io_membership %>%
  filter(year == 1985) %>% 
  ggplot(., aes(formals_sum_joined, informals_sum_joined)) +
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  geom_abline(intercept = 0, slope = 1, color = "grey80") +
  geom_jitter(width = 0.25, height = 0.25, size = 4, alpha = 0.65, shape = 19)+
  labs(title="1985", 
       x = "", y = "") +
  scale_x_continuous(breaks = seq(0, 10, 2), limits = c(0, 10)) +
  scale_y_continuous(breaks = seq(0, 8, 2), limits = c(0, 8)) +
  coord_cartesian(xlim=c(-1,11), ylim=c(-1,11)) +
  coord_fixed(ratio = 1) +
  scatter_theme

set.seed(20210513)
j_1995 <- io_membership %>%
  filter(year == 1995) %>% 
  ggplot(., aes(formals_sum_joined, informals_sum_joined)) +
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  geom_abline(intercept = 0, slope = 1, color = "grey80") +
  geom_jitter(width = 0.25, height = 0.25, size = 4, alpha = 0.65, shape = 19)+
  labs(title="1995", 
       x = "", y = "") +
  scale_x_continuous(breaks = seq(0, 10, 2), limits = c(0, 10)) +
  scale_y_continuous(breaks = seq(0, 8, 2), limits = c(0, 8)) +
  coord_cartesian(xlim=c(-1,11), ylim=c(-1,11)) +
  coord_fixed(ratio = 1) +
  scatter_theme

set.seed(20210513)
j_2005 <- io_membership %>%
  filter(year == 2005) %>% 
  ggplot(., aes(formals_sum_joined, informals_sum_joined)) +
  geom_smooth(fullrange=TRUE, size = 2, color = "grey50") +
  geom_abline(intercept = 0, slope = 1, color = "grey80") +
  geom_jitter(width = 0.25, height = 0.25, size = 4, alpha = 0.65, shape = 19)+
  labs(title="2005", 
       x = "", y = "") +
  scale_x_continuous(breaks = seq(0, 10, 2), limits = c(0, 10)) +
  scale_y_continuous(breaks = seq(0, 8, 2), limits = c(0, 8)) +
  coord_cartesian(xlim=c(-1,11), ylim=c(-1,11)) +
  coord_fixed(ratio = 1) +
  scatter_theme
}

joined_fouryears <- gridExtra::grid.arrange(
  j_1975, j_1985, j_1995, j_2005,
  nrow = 2, ncol = 2,
  bottom=grid::textGrob("Formal IO memberships",
                        gp=grid::gpar(fontsize=30,font=8)),
  left=grid::textGrob("Informal IO memberships", rot = 90,
                      gp=grid::gpar(fontsize=30,font=8)))

ggsave("data/output/figure_app2_joined.pdf",
       joined_fouryears,
       width=16, height=12, units = "in")


#### Appendix standardized coefficients plot ####

## Load in Mansfield and Pevehouse results 

# Tidy
mp_esttab <- read_csv("data/output/mp_esttab.csv", skip = 0L) %>% 
  slice(2:3) %>% 
  select(-1) %>% 
  pivot_longer(names_to = "term", values_to = "value", mp_mp:mp_total) %>% 
  mutate(output = c(rep("beta", 4), rep("se", 4))) %>% 
  pivot_wider(names_from = output, values_from = value) %>% 
  mutate_at(vars(beta:se), ~as.numeric(.)) %>% 
  mutate(conf_low90 = beta - 1.645*se,
         conf_high90 = beta + 1.645*se,
         conf_low95 = beta - 1.96*se,
         conf_high95 = beta + 1.96*se,
         order = seq(4, 1, -1))
mp_esttab

## Load in Greenhill results

# Tidy
greenhill_esttab <- read_csv("data/output/greenhill_esttab.csv", skip = 0L) %>% 
  slice(3:10) %>% 
  select(-1) %>% 
  pivot_longer(names_to = "term", values_to = "value", 
               greenhill_greenhill:greenhill_total) %>% 
  drop_na() %>% 
  mutate(output = rep(c("beta", "se"), times = 4)) %>% 
  pivot_wider(names_from = output, values_from = value) %>% 
  mutate_at(vars(beta:se), ~as.numeric(.)) %>% 
  mutate(conf_low90 = beta - 1.645*se,
         conf_high90 = beta + 1.645*se,
         conf_low95 = beta - 1.96*se,
         conf_high95 = beta + 1.96*se,
         order = seq(4, 1, -1)) 
greenhill_esttab

## Load in Bernauer et al results

# Tidy
bkks_esttab <- read_csv("data/output/bkks_esttab.csv", skip = 0L) %>% 
  slice(3:10) %>% 
  select(-1) %>% 
  pivot_longer(names_to = "term", values_to = "value", 
               bkks_bkks:bkks_total) %>% 
  drop_na() %>% 
  mutate(output = rep(c("beta", "se"), times = 4)) %>% 
  pivot_wider(names_from = output, values_from = value) %>% 
  mutate_at(vars(beta:se), ~as.numeric(.)) %>% 
  mutate(conf_low90 = beta - 1.645*se,
         conf_high90 = beta + 1.645*se,
         conf_low95 = beta - 1.96*se,
         conf_high95 = beta + 1.96*se,
         order = seq(4, 1, -1))
bkks_esttab  

## Combine outputs from different replications

model_combine <- bind_rows(greenhill_esttab, mp_esttab, bkks_esttab) %>% 
  separate(term, into = c("author", "io")) %>% 
  select(-order) %>%
  mutate(rank = ifelse(author == "greenhill", 4, 
                       ifelse(author == "mp", 2, 0))) %>% 
  mutate(rank = ifelse(io == "formal", rank+0.1,
                       ifelse(io == "informal", rank-0.1,
                              ifelse(io == "total", rank-0.3, rank+0.3)))) %>% 
  mutate(io = ifelse(io != "formal" &
                       io != "informal" & 
                       io != "total", "Original", io),
         io = str_to_sentence(io)) %>% 
  mutate(label = ifelse(author == "greenhill", "Greenhill (2010)",
                        ifelse(author == "mp", "Mansfield and\nPevehouse (2006)",
                               "Bernauer\net al. (2010)"))) 
model_combine 

### Plot

model_combine %>% 
  ggplot(., aes(beta, rank, shape = io)) +
  scale_y_continuous(breaks = c(4, 2, 0), 
                     labels = unique(model_combine$label)) +
  labs(x="Standardized regression coefficient\n[90% and 95% confidence intervals]", 
       y="",
       shape = "IO measure") +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill=NA), 
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.minor.x = element_line(colour = "grey90"),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        axis.line = element_line(size = .5, colour = "black")) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black") +
  geom_segment(aes(x = conf_low95, xend = conf_high95, y = rank, yend = rank),
               size = 1, color = "grey80") +
  geom_segment(aes(x = conf_low90, xend = conf_high90, y = rank, yend = rank),
               size = 2, color = "grey50") +
  geom_point(size = 5) +
  scale_shape_manual(breaks=c("Original","Formal","Informal", "Total"),
                     values = c(1,16,18,17))
ggsave("data/output/figure_app3_standardized_coefplot.pdf", units = "in", width = 10, height = 6)

## Use qnorm to find the quantiles where the confidence intervals no longer overlap
model_combine %>% 
  select(author, io, beta, se, conf_low95, conf_high95) %>% 
  filter(io == "Formal" | io == "Informal")

## Define a sequence of values to evaluate
sequence <- seq(0,1,0.0005)

## Mansfield and Pevehouse
#  mean is the beta; sd is the standard error from the regression models
cbind(sequence, 
      formals = qnorm(p = sequence, mean = 0.199, sd = 0.0632, lower.tail  = T),
      informals = qnorm(p = 1-sequence, mean = 0.0340, sd = 0.0498, lower.tail  = T)) %>% 
  as.data.frame() %>% 
  mutate(diff = formals - informals) %>% 
  mutate(min = min(abs(diff))) %>% 
  filter(min == abs(diff))
# 0.072
1-(2*0.072) # the estimates no longer overlap with 85.6% confidence intervals

## Bernauer et al.
cbind(sequence, 
      formals = qnorm(p = sequence, mean = 0.249, sd = 0.0836, lower.tail  = T),
      informals = qnorm(p = 1-sequence, mean = 0.00444, sd = 0.0542, lower.tail  = T)) %>% 
  as.data.frame() %>% 
  mutate(diff = formals - informals) %>% 
  mutate(min = min(abs(diff))) %>% 
  filter(min == abs(diff))
# 0.038
1-(2*0.038) # the estimates no longer overlap with 92.4% confidence intervals




